function [Z, K] = ibphrnd(N, alphav)
H = 2^N-1;
ZZ = zeros(N, H);
m = zeros(1, H);
for h = 1:H
    bcode = dec2bin(h)-'0';
    m(h) = sum(bcode);
    ZZ(N-length(bcode)+1:N, h) = bcode';
end

K = poissrnd((alphav./factorial(N)).*(factorial(N-m).*factorial(m-1)));
Z = zeros(N, sum(K));
kk = 1;
for k = H:-1:1
    if K(k) > 0
        Z(:,kk:kk+K(k)-1) = repmat(ZZ(:,k), 1, K(k));
        kk = kk + K(k);
    end
end

end
